# Code for Adversary_Military and Adversary_Peace Graphs
# Jan 07, 2022

rm(list=ls())
library(ggplot2)
library(Rmisc)
library(dplyr)
library(tidyr)
library(dotwhisker)
library(broom)
library(rtf)
library(readxl)
library(reshape2)
library(haven)
library(naniar)
library(janitor)
library(estimatr)

setwd("~/Dropbox/NikharIndiaSurvey/ReplicationArchive/Dataverse/Data")

us_raw <- read.csv("US_SurveyData_insheet.csv", sep = "")


# Declare Function
interaction_plot_binary <- function(model, effect, moderator, interaction, varcov="default", conf=.95, title="Marginal effects plot", xlabel="Value of moderator", ylabel="Estimated marginal coefficient", factor_labels=c(0,1)){
  
  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = vcov(model)
  }else{
    covMat = varcov
  }
  
  # Extract the data frame of the model
  mod_frame = model.frame(model)
  
  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]
  
  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- c(0,1)
  
  # Compute marginal effects
  delta_1 = beta_1 + beta_3*x_2
  
  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction, interaction] + 2*x_2*covMat[effect, interaction]
  
  # Standard errors
  se_1 = sqrt(var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1
  
  # Determine the bounds of the graphing area
  max_y = max(upper_bound)
  min_y = min(lower_bound)
  
  # Initialize plotting window
  plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(-.5, 1.5), xlab=xlabel, ylab=ylabel, main=title, xaxt="n")
  
  # Plot points of estimated effects
  points(x=x_2, y=delta_1, pch=16)
  
  # Plot lines of confidence intervals
  lines(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), lty=1)
  points(x=c(x_2[1], x_2[1]), y=c(upper_bound[1], lower_bound[1]), pch=c(25,24), bg="black")
  lines(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), lty=1)
  points(x=c(x_2[2], x_2[2]), y=c(upper_bound[2], lower_bound[2]), pch=c(25,24), bg="black")
  
  # Label the axis
  axis(side=1, at=c(0,1), labels=factor_labels)
  
  # Add a dashed horizontal line for zero
  abline(h=0, lty=3)
  
}


require(Matching)
setwd("~/Dropbox/NikharIndiaSurvey/ReplicationArchive/Dataverse/Output")

#Adversary and Military
M3 <- lm_robust(support_trade ~ adversary_treat + military_treat + adversary_treat*military_treat, data=us_raw)
summary(M3)
# Create an interaction plot with binary moderator
pdf(file="Figure3a.pdf", height=6)
adv_mil <- interaction_plot_binary(M3, effect="adversary_treat", moderator="military_treat", interaction="adversary_treat:military_treat", factor_labels=c("No Military Benefit","Military Benefit"), xlabel="Treatment varies whether or not trade will benefit the other country's military", ylabel="Marginal effect of adversary treatment", title="Marginal Effect of Adversary Treatment on Trade Support")

adv_mil
dev.off()
adv_mil



#Adversary and Peace
M5 <- lm(support_trade ~ adversary_treat + peace_treat + adversary_treat*peace_treat, data=us_raw)
summary(M5)
# Create an interaction plot with binary moderator
pdf(file="Figure3b.pdf", height=6)
adv_peace <- interaction_plot_binary(M5, effect="adversary_treat", moderator="peace_treat", interaction="adversary_treat:peace_treat", factor_labels=c("Will Not Help","Will Help"), xlabel="Treatment varies whether or not trade will help ensure peace", ylabel="Marginal effect of adversary treatment", title="Marginal Effect of Adversary Treatment on Trade Support")

adv_peace
dev.off()
adv_peace
